function [X1_save, X2_save, E] = Rulkov(alpha, sigma, mu, k, phi0, N_total, epsilon, N_cut)
% 计算离散忆阻耦合双神经元模型的激发序列及同步误差
%
%   输入参数：
%       alpha, sigma, mu  - Rulkov 映射和模型的参数
%       k                 - 耦合强度
%       phi0              - 忆阻器初始状态
%       N_total           - 总迭代次数
%       epsilon           - phi 的更新系数
%       N_cut             - 去除的瞬态迭代次数
%
%   输出参数：
%       X1_save, X2_save  - 两个神经元的激发序列
%       E                 - 计算得到的同步误差（归一化）
%
% 此函数计算神经元的激发序列，并返回同步误差

    % 初始化神经元状态
    X1 = 0.01;  Y1 = 0;
    X2 = -0.01; Y2 = 0;
    phi = phi0;

    % 初始化保存数据的向量
    X1_save = zeros(1, N_total);
    X2_save = zeros(1, N_total);
    X1_save(1) = X1;
    X2_save(1) = X2;

    % 定义 Rulkov 映射的匿名函数（分段函数）
    rulkovF = @(x, y, alpha) ( (x <= 0) .* (alpha ./ (1 - x) + y) + ...
                               ((x > 0) & (x < alpha + y)) .* (alpha + y) + ...
                               (x >= alpha + y) .* (-1) );

    % 初始化同步误差累加器
    sum_diff = 0;
    N_calc = N_total - N_cut;  % 有效迭代步数

    % 迭代计算
    for n = 1 : (N_total - 1)
        % 计算当前状态下映射输出
        F1 = rulkovF(X1, Y1, alpha);
        F2 = rulkovF(X2, Y2, alpha);

        % 更新神经元状态
        X1_next = F1 + k * (X1 - X2) * tanh(phi);
        Y1_next = Y1 - mu * (X1 + 1) + mu * sigma;

        X2_next = F2 - k * (X1 - X2) * tanh(phi);
        Y2_next = Y2 - mu * (X2 + 1) + mu * sigma;

        % 更新 phi
        phi_next = phi + epsilon * (X1 - X2);

        % 状态赋值
        X1 = X1_next; Y1 = Y1_next;
        X2 = X2_next; Y2 = Y2_next;
        phi = phi_next;

        % 保存当前迭代结果
        X1_save(n+1) = X1;
        X2_save(n+1) = X2;

        % 去除瞬态后的同步误差累加
        if n > N_cut
            sum_diff = sum_diff + abs(X1 - X2);
        end
    end

    % 计算同步误差
    E = sum_diff / N_calc;  % 平均同步误差
end
